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Abstract. 

The dependence of the unfolding pathway of proteins on the pulling speed is 
investigated. This is done by introducing a simple one-dimensional chain comprising N 
units, with different characteristic bistable free energies. These units represent either 
each of the modules in a modular protein or each of the intermediate “unfoldons” 
in a protein domain, which can be either folded or unfolded. The system is pulled 
by applying a force to the last unit of the chain, and the units unravel following a 
preferred sequence. We show that the unfolding sequence strongly depends on the 
pulling velocity Vp. In the simplest situation, there appears a critical pulling speed Vc- 
for pulling speeds Vp < Vc, the weakest unit unfolds first, whereas for Vp > Vc it is the 
pulled unit that unfolds first. By means of a perturbative expansion, we find quite an 
accurate expression for this critical velocity. 



Pulling speed dependence of the unfolding pathway of proteins 

1 . Introduction 


2 


Since the late twentieth century, research on the mechanical stability of macromolecules 
turned a prolihc held due to the advances in manipulation techniques of individual 
biomolecules, usually termed single-molecule experiments. One of the most important 
techniques is atomic force microscopy (AFM), in which a biomolecule is stretched 
between a rigid platform and the tip of the cantilever In these experiments, 

the controlled parameter is either the length of the macromolecule (length-control 
protocols) or the force exerted over it (force-control protocols), and its conjugated 
magnitude is measured. As a result, a force-extension curve (FEC) is obtained, 
which characterizes the elasto-mechanical behaviour of the macromolecule and provides 
fundamental information about its unfolding pathway j5 14 . 


In a typical pulling experiment, the end-to-end distance of the molecule L is 
increased with a pulling rate Vp, that is, dL{t)/dt = Vp. Remarkably, the FEC exhibits 
a sawtooth pattern [9|- p)^ [l^ [Ts] showing how the macromolecule comprises several 
structural units or blocks, in general each one with different stability properties. Each 
block unfolds individually causing a drop in the measured force. The unfolding pathway 
is, basically, the order and the way in which the structural blocks of the macromolecule 
unravel. 

Recent studies show that the pulling velocity plays a relevant role in determining 
the unfolding pathway [14,16-18 . Different unfolding pathways are observed depending 


on (a) which of the ends (C-terminus or N-terminus) the domain is actually pulled from 
and (b) the pulling speed. It has been claimed that it is the inhomogeneity in the 
distribution of the force across the protein, for high pulling speeds, that causes the 
unfolding pathway to change 14,16-18 . Nevertheless, to the best of our knowledge, a 


theory that explains this crossover is still lacking. 

We focus on the unfolding pathway of protein domains composed of several stable 
structural units In this respect, a good candidate is represented by the 

Maltose Binding Protein (MBP), a stable and well characterized protein comprising 
two domains. MBP has been recently employed in studies on mechanical unfolding and 
translocation 17,19-24 . In particular, Bertz and Rief [T^ identihed four intermediate 


states in the FEC of mechanical denaturation experiments, each one associated to 
the unravelling of a specihc unit. The authors termed such units “unfoldons” and 
determined their typical unfolding sequence. However, simulations [^ showed that 
this unfolding scenario holds only for low pulling speed as the pathway depends on the 
rate at which the molecule is pulled: at very low pulling rates, it is the weakest unfoldon 
that unfolds hrst, while at higher rates the hrst unfoldon to unravel is the pulled one. 

The main purpose of this paper is to develop a physical theory that predicts the 
observed dependence of the unfolding pathway on the pulling velocity. Specihcally, we 
do so in the adiabatic limit, that is, the regime of slow pulling speeds that allow the 
system to sweep the whole stationary branches of the FEC 25,27 . Also, we would like 


to stress that, although our approach is applied to the unfolding of a single protein, it 
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extends to the unfolding of polyproteins comprising several domains (modular proteins), 


which unfold sequentially 10 


The plan of the paper is as follows. In section we put forward our model to 
investigate the unfolding pathway of proteins. In general, the free energy characterizing 
each unit (unfoldon or module) of the protein is different, that is, there is a certain 
degree of asymmetry (or disorder) in the free energies. Moreover, we discuss the role of 
thermal noise and the (ir)relevance of the details of the device controlling the length of 
the protein. Section is devoted to the analysis of the pulling of the model, by means 
of a perturbative expansion in both the asymmetry of the free energies and the pulling 
speed. In sections |3.1 and 3.2 we obtain the corrections introduced by the asymmetry 
and the finite value of the pulling speed, respectively. In section 3^ we show that, in 
the simplest situation, there naturally appears a critical pulling speed Vc, below (above) 
which it is the weakest (pulled) unit that unfolds first. In general, when more than one 
unit has a different free energy, we show that for low (high) enough pulling velocity, it 
is still the weakest (pulled) unit that unfolds first, but other pathways are present for 
intermediate velocities. Numerical results for some particular situations are shown in 
section They are compared to the analytical results previously derived, and a quite 
good agreement is found. Finally, section deals with the main conclusions of the paper. 
The appendices cover some technical details that we have omitted in the main text. 


2. The model 


Let us consider a certain protein domain comprising N unfoldons (or a polyprotein 
composed of N, possibly different, modules). From now on, we will refer to these 
unfoldons or modules as units. When the molecule is submitted to an external force 
F, the simplest description is to portray it as a one-dimensional chain, where the end- 
to-end extension of the Tth unit in the direction of the force is denoted by Xj. In a 
real AFM experiment, the molecule is attached as a whole to the AFM device and 
stretched. Following Guardiani et ah |^, we model this system with a sequence of 
nonlinear bonds, as in figure the endpoints of the f-th unit are denoted by gi_i and 
Qi, so that its extension Xi is 

Xi = qi-qi_i, i = l,...,N. (1) 


The evolution of the system follows the coupled overdamped Langevin equations 


IQi = --^G{qi, ..., gjv) + Vi, 
oqi 


( 2 ) 


in which the pi are the Gaussian white noise terms, such that {Vi(t)) = 0 and 
{Vi{t)Vj{'t')) = ‘^ikBTdijdit — t'), with 7 and T being the friction coefficient of each unit 
(the same for all) and the temperature of the fluid in which the protein is immersed, 
respectively {ks is the Boltzmann constant). The global free energy function of the 
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system is 


N 


G{qi ,. . . , ^ Ui{qi - qi-i) + Up{qN) • (3) 

i=l 

In (|^, Up{qN) is the contribution to the free energy introduced by the force-control or 
length-control device (see below), while Ui{xi) is the single unit contribution to G. 

The total length of the system is given by 


N 

^Xi = qN. (4) 

i=l 


In force-control experiments, the applied force F is a given function of time, whereas 
in length-control experiments the device tries to keep the total length g^r equal to the 
desired value L, also a certain function of time. The corresponding contributions to the 
free energy are 

Up{qN) = —FqN, force-control, (5a) 

Up{qN) = - Lf, length-control, (5b) 


in which kp stands for the stiffness of the length-control device. The length is perfectly 
controlled in the limit kp oo, when q^ = L for all times. A sketch of the model is 
presented in hgure 

An apparently similar system, in which each module of the chain follows the 
Langevin equation yxj = —dG/dxi -|- rji has been recently analyzed 27,35 . In this 


approach, the modules are completely independent in force-controlled experiments 
because these Langevin equations completely neglect the spatial structure of the chain. 
While this simplifying assumption poses no problem for the characterization of the force- 
extension curves in 27 , it is not suited for the investigation of the unfolding pathway, in 


which the spatial structure plays an essential role. The spatial structure of biomolecules 
can be described in quite a realistic way by using a model proposed by Hummer and 
Szabo several years ago to investigate their stretching 28 , but the simplihed picture 


which follows from hgure makes an analytical approach feasible. 

Now, we turn to look into the unfolding pathway of this system. As the evolution 
equations are stochastic, this pathway may vary from one trajectory of the dynamics to 
another. Nevertheless, in many experiments 16,17,19 it is observed that there is a quite 
well-dehned pathway, which suggests that thermal huctuations do not play an important 
role in determining it. Physically, this means that the free energy barrier separating the 
unfolded and folded conformations at coexistence (that is, at the critical force, see below) 
is expected to be much larger than the typical energy ksT for thermal fluctuations. 
Therefore, we expect the thermal noise terms in our Langevin equations to be negligible 
and, consequently, they will be dropped in the remainder of the theoretical approach 
developed in the paper. Of course, if the unfolding barrier for a given biomolecule 
were only a few ksTs, the thermal noise terms in the Langevin equations could not be 
neglected and our theoretical approach would have to be changed. 












Pulling speed dependence of the unfolding pathway of proteins 


5 


Uii'Xi) U2{X2) Usixs) U4{x4) Up (q^) 

fA/WVV*AAM/WW\/VW 


Xi Xo X4 Xa 

qo qi q2 qs 


control 
q^ (iovicc 


Figure 1: Sketch of the model for a protein with four units. Each unit is represented 
by a nonlinear spring with potential Ui{xi), in which Xi is the unit’s extension. The 
beads mark the coordinates of their endpoints, so that Xi = — gj_i (by dehnition, 

qo = 0). Finally, the rectangle stands for the device attached to the pulled end q^, 
which controls either the force applied to the molecule (force-control) or its end-to-end 
distance (length-control). The contribution of this device to the system’s free energy is 
Up{qA)-i as shown by ([^ and (|^. 


In order to undertake a theoretical analysis of the stretching dynamics, one further 
simplihcation of the problem will be introduced. We consider that the device controlling 
the length is perfectly stiff, and the total length qj^ = L does not fluctuate for all times. 
We expect this assumption to have little impact on the unfolding pathway: otherwise, 
the latter would be more a property of the length-control device than of the chain. In 
fact, we show in section]^ that the unfolding order is not affected by this simplihcation. 
For perfect length-control, the mathematical problem is identical to that of the force- 
control situation, but now the force F is an unknown (Lagrange multiplier) that must 
be calculated at the end by imposing the constraint qN = Therefore, the 


extensions xfs obey the deterministic equations 

jXi = -U[{xi) + U{X2), (6a) 

yij = -2U[{xi) + U[j^^{xi+i) + U[_^{xi_i), I <i<N, (6b) 

yxAT = -2U'^{xn) + U'^_i{xN-i) + E, (6c) 

F =-iVp + U'j^{xN)- (6d) 

We have introduced the pulling speed 

Vp = L, (7) 


which is usually time independent. 

We assume that Ui{xi) allows for bistability in a certain range of the external force 
F, in the sense that Ui{xi) —Fxi is a double-well potential with two minima, see hgure 
1^ Therefore, in that force range, each unit may be either folded, if Xi is in the well 
corresponding to the minimum with the smallest extension, or unfolded, when Xi belongs 
to the well with the largest extension. If the length is kept constant [vp = 0), there is 
an equilibrium solution of ([^, 

U[{Ff) = U'ixf) = ■■■ = U'^{x%) = F^\ 


( 8 ) 
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and F®* is calculated with the constraint This solution is stable as long as 

U”{xf) > 0 for all i. 

If all the units are identical, Ui{x) = U{x), the metastability regions of each module 
(the range of forces for which the equation U[{x) = F has several solutions) coincide. 
Therefore, we obtain stationary branches corresponding to J unfolded units and N — J 
folded units that have been analyzed in detail in 25,27. If all the modules are not 


identical, the metastability regions do not perfectly overlap since the units are not 
equally strong: the weakest one is that for which the equation U[{x) = F ceases to have 
multiple solutions for a smaller force value. 

It is important to note that we can change all the forces U[{xi) to Vl^xf) = 
Ul{xi) — Fq and F to cp = F — Fq, we have the same system (§ but with Vf and 
(p instead of U'^ and F, respectively. Then, we may use the free energies for any common 
value of the force Fq and interpret the Lagrange multiplier as the excess force from this 
value to be applied to the system[U 


3. The pulled chain 

Let us consider the pulling of our system. We write the Tth-unit free energy as 

Ui{x) = U{x)+f5Ui{x), (9) 

in which U{x) is the “main” part, common to all the units, and f6Ui{x) represents 
the separation from this main contribution. If all the units are perfectly identical, 
Ui{x) = U{x) for all i or, equivalently, 6Ui{x) = 0. In principle, in an actual experiment, 
the splitting of the free energy in ([^ can be done if the free energy Ui of each unit is 
known: we may dehne the common part as the “average” free energy over all the units, 
U{x) = U{x) = N~^ Ui{x), and f5Ui{x) = Ui{x) — U{x). From a physical point of 
view, the dimensionless parameter ^ > 0 measures the importance of the heterogeneity 
in the free energies. Our theory could be applied to a situation in which the free 
energy deviations 5Ui were stochastic and followed a certain probability distribution. 


for instance to represent the slight differences among very similar units, as done in 27 


to analyze the force-extension curves. In particular, the forces U[{x) in the evolution 
equations can also be split as 

U[{x) = U\x) + f 6fi{x), 6fi{x) = 6U-{x). (10) 

As already noted above, we can use the free energies for any common value of the 
force Fq, and interpret F as the extra applied force from this value. In what follows, 
we consider U{x) with two, equally deep, minima corresponding to the folded (F) and 
unfolded (U) conhgurations. Figure presents a qualitative picture of the free energy 
and its derivative. The two minima correspond to lengths £f and iu, with ip < ip. 
Also the point 4 at which U”{ib) = 0 is marked. 

f A similar result is also found if the length is controlled by using a device with a finite value of the 
stiffness kp. A constant force only shifts the equilibrium point of a harmonic oscillator: (g^v ~ L) must 
be substituted by (g^r — L — Fo/kp). 
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Figure 2: Qualitative behaviour of the main contribution to the free energy U{x) at the 
critical force (top panel) and associated force U'{x) (bottom panel) as a function of the 
extension. The values of the lengths at the folded and unfolded minima are ip and ip, 
respectively, whereas the threshold length ij, stands for the length corresponding to the 
limit of stability and F), is the corresponding force. 


It is the condition U"{ih) = 0 that essentially determines the stability threshold, 
as it provides the limit force F), = U'{i}f) > 0 at which the folded basin ceases to exist 
for the “main” potential. In the deterministic approximation considered here, thermal 
fluctuations are neglected and, for F < F^, the folded unit cannot jump over the free 
energy barrier hindering its unfolding: it has to wait until, aX F = F),, the only possible 
extension is that of the unfolded basin. Of course, neglecting thermal noise restricts in 


some way the range of applicability of our results, see section 3.3 for a more detailed 
discussion and also the numerical section HI 

Keeping the above discussion in mind, now we analyze the limit of stability of 
the different units. The asymmetry correction 6fi shifts the threshold force for the 
different units. The extension Xi^b at which the Fth unit loses its stability is obtained 
by solving the equation U”{xi^b) = U"{xi^b) = 0, which linearized in both the 

displacement Xi^b — h and f reads 


U"{ib) + U'"{ib){x,,b - h) + ^df'ih) 
Noting that U"{ib) = 0, we get that 

SfliQ 


= 0 . 


( 11 ) 


^i,b ib ^ 


U"'{ib 


( 12 ) 
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See [Appendix A for details. The corresponding force 


IS 


F,,t = F, + ^6fM, 


( 13 ) 


in which we have also dropped terms of the order of Then, units with 6fi{if,) < 0 
{6fi{ib) > 0 ) are weaker (stronger) than average. 

When the system is continuously pulled, the total length of the system L has 
been shown to be a good reaction coordinate [^. Therefore, on physical grounds it is 
reasonable to use L to measure time and write the evolution equations (|^ as 

(jx-i 

IVp— = -U'{xi) + U\x 2 ) + ^[-5fi{xi) + Sf 2 {x 2 )], (14a) 

dx ■ 

= -‘^U'{xi) + U'{xi+i) + U'{xi_i) + f[-26fi{xi) + 6fi+i{xi+i) + 6fi_i{xi_i)], 

1 <i < N, (14b) 

= -2U\xn) + U'{xn-i) + F + ^[-26fN{xN) + 6fN-i{xN-i)], (14c) 

F = 'yVp + U'{xN) +^dfN{xN)- (14d) 


Moreover, this change of variable makes the pulling speed Vp appear explicitly in the 
equations, allowing us to consider Vp as a perturbation parameter for slow enough pulling 
processes. 

Now, we consider a system in which the asymmetry in the free energies is small 
and which is slowly pulled. Thus, (14) is solved by means of a perturbative expansion 
in powers of the pulling velocity Vp and the disorder parameter that is. 


Xi{L) = xf\L) + ^6xi{L) + VpAxi{L), (15a) 

F{L) = F(°)(L) + ^6F{L) + VpAF{L), (15b) 

up to the linear order in both Vp and f. 

The zero-th (lowest) order corresponds to the chain of identical units (^ = 0) with 
a given constant length L {vp = 0). Namely, xf'^ and F^^^ obey the equations 

t) = -U'{xf) + U'{xf), (16a) 

0 = -2U\xf'^) + U\xf}f) + U\xf\), 1 < * < AT, (16b) 

0 = -2U\x%'>) + f/'(a;5;Li) + F^^\ (16c) 

FW = f/'(a:®). (16d) 

The solution of this system is straightforward, 

U'ixf^) = (17) 


the force is equally distributed among all the units of the chain in equilibrium, as 
expected. If we start the pulling process from a conhguration in which all the units are 
folded and the force is outside the metastability region (the usual situation), the units 
extensions and the applied force are 



L 

iv’ 


V'i, F<“l = (/'(<), 


( 18 ) 
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to the lowest order. To calculate the linear corrections in ^ and Vp, we have to substitute 


(15) and (18) into (14), and equate terms proportional to f and Vp, respectively. This 


is done below in two separate sections: firstly, for the asymmetry contribution 5xi and, 
secondly, for the “kinetic” contribution Ax*. 


3.1. Asymmetry correction 

All the modules are not characterized by the same free energy, and here we calculate 
the first order correction introduced thereby. The asymmetry corrections 6xi obey the 
system of equations 

5fi{^) - 


6x2 — 6x1 = 


U"{i) 


Sxi+i + 6xi_i — 26xi = 


2dMe) - 5h 


i+l I 




SxN-^ — 2 Sxn = 


U"{i) 

26fNii) - 6fN-iii) - 6F 
U"{i) 


(19a) 

l<i<N (19b) 

(19c) 

(19d) 


6F = U'\i)6xN + SfN{i), 

which is linear in the 6xfs, and thus can be analytically solved. It is clear that our 
expansion breaks down when U”{£) = 0. This was to be expected, since we know that 
the stationary branch with all the modules folded is unstable when Uf becomes negative 
for some unit i, and to the lowest order this takes place when U"{€) = 0. 

The solution of the above system of difference equations is obtained by standard 
methods (sTl, with the result 


N 


= vj , sF = sm = ifY.^m. 


u"{e) 


( 20 ) 


2 = 1 


See Appendix B for more details. Interestingly, the force is homogeneous across the 
chain, since to first order in ^ we have that 

f/'(a;,) = U\xf^) + f[U^\xf)Sx, + d^xf)] = U'{£) + (Ef{£) = + ^dF (21) 


This is nothing but the stationary solution ([^, up to first order in the disordeij§j 
Moreover, (20) implies that there are units with dxi > 0 and others with dxi < 0, 


depending on the sign of df{£) — dfi{i). This is a consequence of our perturbation 
expansion, since = L for all times, as given by (j^, and thus dxi = 0. 

Let us remember that we denote by the value of the extension at which the 
common main free energy reaches its limit of stability, see figure Taking into account 
only the asymmetry correction, it is the weakest unit that unfolds first, since the most 
negative dfi{i) leads to the largest positive dxi and then it is the one that first verifies 

§ If the zero-th order free energy were the average of the th’s, no correction for the Lagrange multiplier 
(applied force) would appear to the first order. This is logical, up to the first order the force expression 
coincides with the spatial derivative of the average potential, that is, fF)-\-^SF = Ufi)+£Sf{£) = 
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the condition Xi = i ^Sxi = ib (for a more detailed discussion, see Appendix A). An 


alternative way of looking at this is to recall that the force corresponding to the limit of 
stability is smallest for the weakest unit: since the force is homogeneously distributed 
along the chain, it is the weakest module that hrst reaches its stability threshold. 


3.2. Correction due to the finite pulling speed 

Now we look into the “kinetic” correction due to the hnite pulling speed Vp. The zero-th 
order solution is given by (18), so that dxf'^/dL = N~^ for all i, and we have 


Axo — Axi = 


7 


NU"{iy 

Axj+i + Axi-i — 2Axi = 

1 


7 


AxAr_i — 2Axn = 

AF = 7 




1 <i < N, 


U"{tj 

U"{i)AxN. 

The solution to this system of linear difference equations (Ml is 


Axi = 


AF = 


7 


2NU"{i) ^ 

(Ar + l)(2Ar + l)7 

6N 


{N + 1){N-1) 


(22a) 

(22b) 

(22c) 

(22d) 

(23a) 

(23b) 


See Appendix B for more details. Again, ^2- Axt = 0 because the zero-th order solution 
(18) gives the total length, times. (23a) is reasonable on intuitive 


grounds: the kinetic correction Axi increases with i because the last module is the one 
that is actually pulled. Therefore, on the basis of only the kinetic correction, it is the 
last module that would unfold hrst because Ax^ is the largest. Thus, the condition 
Xi = i + VpAxi = ib is hrst verihed for i = N. 

It is interesting to note that the force was equally distributed for the asymmetry 


correction, as expressed by (21), but this is no longer true if we incorporate the kinetic 
correction. Up to the the hrst order, U[{xi) = U'{xf'’ + f6xi + VpAxf) + ^Sfi{xi) ~ 
U'{i) -|- ^Sf{i) + VpU”{i)Axi. Therefore, the force U[{xi) depends on the unit v. for all 
times, it is smaller the further from the pulled unit we are. Again, there is an alternative 
way of understanding why the last unit would unfold hrst if we were considering perfectly 
identical units (.^ = 0): for any time, it would be the last unit that suhered the largest 
force and thus the hrst that reached their common limit of stability Fb. 


3.3. The critical velocities 

If the last unit is not the weakest, there is a competition between the asymmetry and 
the kinetic corrections. For very low pulling speeds, in the sense that Vp/f —)■ 0, the 
term proportional to Vp can be neglected and it is the weakest unit (the one with the 


largest 6xi) that unfolds hrst, as discussed in section 3.1 On the other hand, for very 
small disorder, in the sense that ^/vp —)■ 0, the term proportional to f is the one to be 



















Pulling speed dependence of the unfolding pathway of proteins 


11 


neglected and it is the last unit (the one with the largest Axi) that unfolds hrst, as also 
discussed in section |3.2[ Therefore, different unfolding pathways are expected as the 


pulling speed changes. 

Collecting all the contributions to the extensions, we have that 


Xi = i + 


f6f{e) - Vp-f 


- 1 
6iV 


U"{i) 


+ 




U"{i) 


(24) 


We have rearranged the terms in x* in such a way that the first two terms on the rhs are 
independent of the unit i, all the dependence of the length of the module on its position 
across the chain has been included in the last term. Note that we are expanding the 
solution in powers of Vp around the “static” solution, which is obtained by putting Vp = 0 


in (24). Thus, the “static” solution corresponds to the stationary one the system would 


reach if we kept the total length constant and equal to its instantaneous value at the 


considered time. It is essential to realise that (24) is only valid for very slow pulling. 


as long as the corrections to the “static” solution are small, and this is the reason why 
the limit of stability is basically unchanged as compared to the static case. In order to 
be more precise, we refer to this kind of very slow pulling as adiabatic pulling. A main 
result of our paper is that, even for the case of adiabatic pulling, there appear different 
unfolding pathways depending on the value of the pulling speed. 

In the adiabatic limit we are considering here, the pulling process has to be slow 
enough to make the system move very close to the stationary force-length branches, 
but not so slow that the system has enough time to escape from the folded basin. 
As discussed in 27 , there is an interplay between the pulling velocity and thermal 


fluctuations. For very slow pulling velocities, the system has enough time to surpass 
the energy barrier separating the two minima, which leads to the typical logarithmic 
dependence of the “unfolding force” Fu on the pulling speed, specihcally Fjj oc (lnxp)“ 
mii On the other hand, as already argued at the beginning of section]^ for adiabatic 


pulling, the units unfold not because they are able to surpass the free energy barrier but 
because the folded state ceases to exist at the force Ff, corresponding to the upper limit 
of the metastability region. 

The unit that unfolds hrst is the one for which x* = £b for the shortest time. In 
light of the above, it is natural to investigate whether it is possible to determine which 
module is the hrst to unfold for a given pulling speed. To put it another way, we would 
like to calculate the “critical” velocities which separate the velocity intervals in which a 
specihc module unfolds hrst. Let us assume that, for a given pulling speed Vp, it is the 
Ath module that unfolds hrst. All the modules j to its left, that is, with j < i, will not 
open hrst if the pulling velocity is further increased because the diherence between the 
kinetic corrections Axj — Axj increases with Vp. Therefore, the hrst module j to unfold 
when the velocity is sufficiently increased it is always to its right. The velocity x*(j) for 


The parameter a is of the order of unity, its particular value depends on the specific shape of the 


potential (linear-cubic, cuspid-like, ...) considered 29 . 
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which each couple of modules (i, j), j > i, reach simultaneously the stability threshold 
is determined by the condition 




(25) 


(25) determines both the value of ic (or time tc) at which the stability threshold is 
reached and the relationship between Vp and f. (24) implies that 


- (SMQ + = -(SMC) + 


2N 


2N 


(26) 


We already know that the length corresponding to the limit of stability is very close 
to the threshold length i/,, its distance thereto being of the order of a/^, as shown in 


Appendix A[ Therefore, to the lowest order, 4 can be approximated by £b, and we get 

7u*(j) _ 2N[6f,{ib) - 6fM] 




jU - 1 ) -*(*- 1 ) ’ 


J >1. 


(27) 


Clearly, the minimum of these velocities is the one that matters: Let us denote by j 
the position of the module for which u*(j) reaches its minimum value uJnin, 


■W 

min 


= = ramv\j) 


(28) 


for Vp just below it is the Ath module that unfolds first, but for Vp just above 
it is the j-th module that unfolds first. Let us denote the weakest module by ai, that 
is, Sfi{ib) is smallest for i = ai. If Up is smaller than the first unit to reach the 
stability limit is the weakest one. Then, we rename the latter velocity that is. 


(5/ai(4) = min(5/i(4), (29) 

I 

because it is the first one of a (possible) series of critical velocities separating different 
unfolding pathways, see below. 

Let us denote by 0)2 the module which unfolds first in the “second” velocity region, 
Vp just above vi^\ that is, a 2 = ■ This unit ceases to be the first to unfold for the 

velocity 

= «3n (30) 

The successive changes on the unfolding pathway take place at the critical velocities 


C ’ = Cn. (31) 

in which This succession ends when ak+i = N: in that case, for Vp > v^\ 

the first unit to unfold is always the pulled one. This upper critical velocity can be 
computed in a more direct way. 


V 


end 

c 


= max'i;A„(Ar). 

3 


(32) 
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Consistency of the theory reqnires that a short proof is presented in 

Appendix C 


We have a trivial case for ai = N, when the pnlled nnit is precisely the weakest and 
it is always the hrst to nnfold for any pnlling speed. The simplest nontrivial case appears 
when all the modnles has the same free-energy with the exception of the weakest, and 


Q!i 7 ^ N, (27), (29) and (32) rednce to 


'yv, 


( 1 ) 


^yend 






2N[6fM{h) - 5fa,{h)] 
N{N — 1) — ai(ai — 1) 


(33) 


Note that the sitnation is qnite simple, since there exist a single critical velocity 


(1) 6iid 

Vc = vy = n® . 


For Vp < Vc 


the weakest modnle nnfolds hrst whereas for Vp > Vc 


the last one unfolds hrst. In general, when the units have diherent free energies the 
situation may be more complex, as shown in the previous paragraph. There appear 
intermediate critical velocities, which dehne pulling speed windows where neither the 
weakest unit nor the last one is the hrst to unfold. In order to obtain these regions, we 


need to recursively evaluate (31). 


4. Numerical results 


Throughout this numerical section, we check the agreement between our theory and 
the numerical integration of the evolution equations. Firstly, we discuss the validity of 
the simplihcations introduced in the development of the theory, namely (i) negligible 
thermal noise and (ii) perfect length-control. Secondly, we look into the critical pulling 
speed, showing that there appears such a critical speed in the simulations and comparing 
this numerical value to the theory developed before. 

We consider a system composed of = 4 unfoldons, such as the maltose binding 
each one characterized by a quartic bistable free energy. In reduced 


17 


protein 

variables, the free energies have the form Ui{x) = eiU{x), where 

U(x) = 1 [[x - a'f - , 


(34) 


with Cj = 1 for i 7 ^ 1, Cl < 1, cr = 0 and a = €D The value of the friction coefficient is, 
also in reduced variables, 7 = 1 . We use these dimensionless reduced variables in order 


to make it easier to compare our results to those in 17 . The function (34) is one of the 


simplest, but reasonable, choice to describe the free energy of different unfoldons of the 
same protein domain. Using our notation, we have 


Sfi(x) = 0 , ! 7 ^ 1 , i5/i(x) = -iu'(x), 


(35) 


^ Here, the value of a is different from the one in (tr = 8). Its only effect is a shift of the origin 
of the extensions, our choice implies that a positive (negative) sign of the extension corresponds to an 
unfolded (folded) configuration. 
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with .^ = 1 — ei- 
asymmetry 


(12) and (13) give us the limits of stability up to first order in the 


^-Fb, i ^ 1, Fi,b ^ (1 - ^i)Fb 


(36) 


For this simple example, (36) is exact. The weakest unit is the first one, because Fi ;, is 
the minimum value of the force at the limit of stability. For the values of the parameters 
we are using, ii, = aj a/ 3 = 1.73 and Fh = U'{ih) = 2\/3a^/9 = 10.4. Since we are writing 
the free energies for a common given value of the force, we are assuming that all the 
units have their two minima equally deep at the same force. This assumption is made to 
keep things simple: the main ingredient for having an unfolding pathway that depends 
on the pulling speed is to have different values of the forces at the stability threshold 
for the different units. 

In the case we are considering, the weakest unit is the first one, while the others 
share the same free energy. This means that we have the simplest scenario for the 
critical velocity in our theoretical approach: it is always the weakest (for Vp < Vc) or the 
last (for Vp > Vc) unit that opens first, as discussed at the end of the previous Section. 


Here, (33) for ai = 1 and iV = 4 reduces to 


7^c _ 2 

T ■ 


(37) 


To start with, we consider the relevance of the noise terms in ([^. In figure 
we plot the integration of the Langevin equations together with the deterministic 
approximation for a concrete case: the free energy of the first unit corresponds 
to ei = 0.8 (.^ = 0.2), the stiffness of the device controlling the length is kp = 5, 
the temperature is T = 1 and the pulling speed is Vp = 0.38. For these values of 


the parameters, taken from 17 , the critical velocity in (37) is Vc = 1.4, so we are 


considering a subcritical velocity, Vp < Vc- Thermal fluctuations are small, and thus 
the same unfolding pathway is observed in the deterministic and the majority of the 
stochastic trajectories. 

Let us consider in more detail the relevance of thermal noise: from a physical point 
of view, it may be inferred by looking at the height of the free energy barrier at the 
critical force in terms of the thermal energy ksT. For the values of the parameters we 
are using, this barrier is around 20 in reduced units. This explains why thermal noise 
is basically negligible in figure in which T = 1. If the temperature is decreased to 
T = 0.25, the barrier is so high, around 80 times the thermal energy, that essentially all 
the stochastic trajectories coincide with the deterministic one. On the other hand, if the 
temperature is increased to T = 4, the barrier in only a few, around 5, times the thermal 
energy, and we expect that the deterministic approximation ceases to be valid. In order 
to further clarify this point, we present figure]^ Both panels display bar graphs with the 
frequencies with which each unit unfolds first in the stochastic trajectories obtained over 
1000 trajectories of the Langevin equations (|^ with perfect length control and different 
values of the temperature. In the left panel, a subcritical velocity Vp = 0.38 < Vc is 
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Figure 3: Evolution of the extensions of the different units in a pulling experiment as a 
function of the length of the system q^. The pulling speed is Vp = 0.38 and the length- 
control device has a stiffness kp = 5. The symbols correspond to a typical realization of 
the Langevin process ([^ with T = 1, whereas the lines correspond to the deterministic 
(zero noise) approximation. 




Figure 4: Frequency with which each of the units unfolds hrst when the Langevin 
equations with perfect length-control are integrated for different values of the 
temperature. (Left) Numerical frequencies obtained in 1000 trajectories, for a subcritical 
pulling speed Vp = 0.38 < Vc- (Right) The same as in the left panel, but for a 
supercritical pulling speed Vp = 2 > Vc- As the temperature decreases, the frequency of 
the deterministic unfolding pathway approaches unity in both cases. 
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Figure 5: (Left) Evolution of the extensions of the different units in a pulling experiment 
as a function of the length of the system q^. The symbols correspond to the integration 
of the deterministic equations, for kp = 5 (hlled symbols) and kp = 50 (empty symbols), 
whereas the line correspond to the limit kp ^ oo. The pulling speed is the same as in 
£gure[^ that is, Vp = 0.38. (Right) Comparison between the desired and actual lengths, 
L and qn, for the different values of the stiffness considered in the top panel. It is 
observed that the length control improves as kp increases. 


considered, so that the weakest (hrst) unit is expected to unfold hrst. In the right 
panel, the numerical data for a supercritical velocity Vp = 2 > Vc are shown, for which 
the pulled (fourth) unit would unfold hrst. The effect of thermal noise is quite similar 
in both cases. For the low temperature T = 0.25, the frequency of the deterministic 
pathway is close to unity and, for the temperature in hgure T = 1, its frequency is 
still very large, clearly larger than any of the others. On the other hand, for the higher 
value of the temperature, T = 4, thermal noise can no longer be neglected. 

In the following, we restrict the analysis to the physically relevant case in which 
the deterministic approximation gives a good description of the hrst unfolding event. 
In hgure(left panel), we look into the same pulling experiment as before, but now we 
compare the deterministic evolution of the extensions for two hnite values of the stihness 
to the kp ^ oo limit. Consistently with our expectations, the unfolding pathway is not 
ahected by this simplihcation. Nevertheless, the control of the length of course improves 
as kp increases (see right panel). Although for the smaller values of kp the length is not 
perfectly controlled, the curves in the left panel, which correspond to diherent values of 
kp, are almost perfectly superimposed when plotted as a function of the real length of 
the system (but not of the desired length L). This means that the real length g^r is 
a good reaction coordinate, as already said in section 

We have integrated the deterministic approximation ([^ (zero noise) of the Langevin 
equations for diherent values of the pulling speed, and extracted from them the numerical 
value of the critical velocity as a function of the asymmetry ^ = 1 — ei. In order to 
obtain this numerical prediction, we initially set Vp equal to the theoretical critical 
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Figure 6: Phase diagram for the unfolding pathway in the pulling velocity-asymmetry 
plane. Two well-dehned zones are separated by the curve giving the critical velocity Vc as 
a function of the asymmetry ^ in the free energy of the first unit. The numerical values 
for Vc (circles) are compared to the theoretical expression (solid line), ( |3^ . The dashed 
line corresponds to the alternative approach discussed in the text, which improves the 
agreement between the numerical results and the theory for ^ > 0.1. Error bars have 
been omitted because they are smaller than point size. 


velocity given by (37). Then, we recursively shift it by a small amount 6vp, such that 
dvp/vc = 0.0001, until the pathway changes. We compare the values so obtained to the 
theoretical expression (37), in figure]^ We find an excellent agreement for ^ < 0.1, for 
^ > 0.1 there appears some qnantitative discrepancies. These discrepancies stem from 
two points: (i) the perturbative expansion used for obtaining (33) from (25) and (ii) the 
intrinsically approximate character of (25), since ib gives rigorously the limit of stability 
only for the static case Vp = 0. Therefore, we have looked for the solution of (25) in the 
nnmerical integration of the deterministic equations. This is the dashed line in figure 
which substantially improves the agreement between theory and nnmerics because we 
have eliminated the deviations arising from point (i) above. In fact, for the case we have 
studied in the previous hgures, which corresponds to a not so small asymmetry ^ = 0.2, 
the improved theory gives an almost perfect prediction for the critical velocity. 

Now we consider a more realistic potential for the nnits, which has been introduced 
by Berkovich et al. for modelling the nnfolding of single-nnit 127 and ubiquitin proteins 
in AFM experiments 33,34 . Moreover, it has also been used to investigate the stepwise 


unfolding of polyproteins in force-clamp conditions 35 and their force-extension curves 
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U{x) = f/o [(1 - ^ 


I 1 X 2x^ . 

+ + ( 38 ) 


that is, it is the sum of a Morse and a worm-like-chain potentials, representing the 
enthalpic and the entropic contributions to the free energy, respectively |33l|34] . We take 


the values of the parameters from 27,33 , P = 0.4nm (persistence length), Lc = 30nm 
(contour length), T = 300K, Uq = lOOpN nm(~24fcsT), = 4nm, 6 = 2. We 
measure force and extensions in the units [F] = 100 pN, Lc = 30 nm, respectively. 
Accordingly, dimensionless variables are introduced with the dehnitions /i = Uo/{Lc[F]), 
/3 = 2bLc/Rc, p = Rc/Lc, A = ksTLc/{4:PUo). Thus, a dimensionless potential is 
obtained, which reads 


f/(x) = /i| [l — e ^ — 1 -t- A -1 — X -1- 2x^^ I , 


(39) 


Note that, in order not to clutter our formulae, we have not introduced a different 
notation for the dimensionless potential. The values of the parameters therein are 
p = 0.0333, /3 = 30, p = 0.133 and A = 0.776. In dimensionless variables, Ff, = 0.527 
(52.7pN) and ih = 0.157 (4.70nm). The relevant time scale is set by the friction 
coefficient 7, [f] = 'yLc/[F]. In turn, 7 is given by the Einstein relation D = ksTh, 
where D is the diffusion coefficient for tethered proteins in solution. We consider a 


typical value of D, also taken from 33 , D = 1500 nm^/s, so that 7 = 0.0028pN nm“^ 

s. 

We consider a system of 4 units, again with all the units but the hrst being identical. 
Then, Ui{x) = U{x), i ^ 1, and the hrst unit is the weakest because Ui{x) = (1—.^)[/(x). 
The situation is then similar to the one we have already analyzed with the quartic 


potential (34), but there is a difference that should be noted: here, Lf{x) is the free 
energy at zero force, whereas for the quartic potential Lf{x) was the free energy at a 
force Fq for which the folded and unfolded minima were equally deep. Then, the force 
here must not be interpreted as the extra force from Fq, but as the whole force that is 
applied to the polyprotein. On the basis of our theory, we expect the simplest situation 
with only one critical velocity Vc, below (above) which the weakest unit (the pulled unit) 
unfolds hrst. This is also indeed the case in the numerical simulations, and we compare 
the theoretical and numerical critical velocities in hgure A very good agreement is 
found again, up to values of the asymmetry ^ of the order 0.1 — 0.2. 

The above discussion shows that the validity of the theory presented here is not 
restricted to simple potentials like the quartic one; on the contrary, it can be conhdently 
applied to experiments in which the units are described by realistic potentials. Moreover, 
for the typical parameters we are using, the theoretical critical velocity Vc for the 
Berkovich potential equals 1270nm/s for an asymmetry ^ = 0.1, which can be regarded 
as quite a conservative estimate of the largest asymmetries for which our theory gives 
an almost perfect description of the unfolding pathway. Interestingly, this pulling speed 
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Figure 7: Phase diagram for the unfolding pathway in the pulling velocity-asymmetry 
plane for the Berkovich potential, (39). Again, as in hgure there appear two well- 
defined pulling regimes, separated by the curve giving the critical velocity Vc as a function 
of the asymmetry The numerical values for Vc (circles) compare very well with the 
theoretical expression (solid line), (33). Again, the dashed line corresponds to the 
alternative approach discussed in the text for the quartic potential, which once more 
signihcantly improves the agreement theory-simulation for the larger values of f. 


corresponds to the upper range of velocities usually employed in AFM experiments, for 
instance see Table I of [^. Therefore, testing our theory in real AFM experiments with 
modular proteins should be achievable. 

Finally, we consider a more complex situation, in which more than one unit is 
different from the rest and there may exist more than one critical velocity. To be 
concrete, we have considered a system with 4 units in which U 2 {x) = Uz{x) = U{x), 
Ui{x) = {l — f)U{x) as before but Ui{x) = (1-|-3.^/2)f/(x), with 17(a;) being the quartic 
potential in (34). In this situation, we have two different critical velocities: for very low 
pulling speeds, the weakest unit is the first to unfold, but there appears a velocity window 
inside which neither the weakest nor the pulled unit is the first to unfold. This stems 
from the fact that the hrst and the third unit reach simultaneously the limit of stability 
for a velocity n^(3) = 4,^7“^Fb/3 that is smaller than the velocity n^(4) = 5f'y~^Fh/3 
for which the first and the last would do so. The physical reason behind this is the 
threshold force of the pulled unit being larger enough than that of the third one. We 
recall that n*(j) is the velocity for which the Ath and the j-th unit reach simultaneously 
their limits of stability. Afterwards, the third unit and the fourth attain the limit of 
stability in unison for a velocity n^(4) = and the following picture emerges 
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from our theory. Using the notation introduced in section |3.3[ we dehne two critical 
velocities, 








= 2K, 


such that for Vp < Vc^\ it is the weakest unit that unfolds hrst, for 

,( 2 ) 


<Vp<Vc 


( 2 ) 


(40) 


it is 


the third unit that unfolds hrst, and, hnally, for Vp > Vc , the hrst unit to unfold is the 
pulled one. 

We check the more complex scenario described in the previous paragraph in hgure 
^ Therein, we observe that (i) our theory correctly predicts the existence of the three 
pulling regimes described above but (ii) even for very small asymmetries, there appear 
some noticeable discrepancy between theory and simulation. Since the validity of the 


perturbative expansion for obtaining the critical velocities from the condition (25) is 


strongly supported by the accurateness of the theoretical prediction for the simplest 
case, see hgures[^and[^ this discrepancy should stem from the intrinsically approximate 
character of the condition U-' = 0 for determining the stability threshold when Vp ^ 0. 
Therefore, an improvement of the present theory should involve the derivation of a 
more accurate condition for obtaining the stability threshold in the case of hnite pulling 
velocity. This point, which probably makes a multiple scale analysis necessary for lengths 
close to the condition t/" = 0, certainly deserves further investigation. 

5. Conclusions 


We have investigated the general properties of the unfolding pathway in pulled proteins 
by means of a simple model portraying the chain as a sequence of nonlinear modules. The 
deterministic approximation of the Langevin equations controlling the time evolution of 
the units extensions is used therein. This is a sensible approach in our model system: 
note that the free energies characterising the different units are considered to be very 
similar, and therefore their Kramers rates for thermal activated unfolding would also be 
very close. Then, were thermal effects important, our unfolding trajectories would be 
essentially stochastic and we would not observe a specihc unfolding pathway. 

Nevertheless, in some recent optical tweezers experiments in which thermal 
fluctuations are relevant, dehnite pathways have also been observed. Although 


trajectories are indeed stochastic, the pathway is well-characterised 45,46 . Therefore, 


the analysis of these experiments needs a more sophisticated theory, which takes into 
account thermal noise effects. Moreover, it seems that there are other elements that 
should be incorporated, such as (i) the possible coupling between the different units and 
(ii) their dissimilar free energies. For example, the former may explain the existence of 


dead ends observed in 46 , that is, intermediate states that do not allow the system to 


completely relax, whereas the latter may lead to Kramers rates leading to the separation 
of the timescales for the different unfolding events. 

The equilibrium extensions of the units are governed by different free energies, 
which we have called asymmetry or disorder in the free energies. Our theory, based 
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Figure 8: Phase diagram for the unfolding pathway in the pulling velocity-asymmetry 
plane for the more complex case discussed in the text. Now, we have three well-dehned 
pulling regimes, separated by two curves giving the critical velocities and as a 
function of the asymmetry In this case, our theoretical approach is able to reproduce 
the existence of the three different pulling regimes, but the discrepancies between the 
theoretical and the numerical values for the critical velocities are larger than in the 
simpler cases considered in the previous figures. 


on stability considerations, is able to explain the experimental observations: (i) for low 
pulling speeds, it is the weakest unit that unfolds hrst, (ii) for large enough pulling speed, 
it is the pulled unit that opens first. This has been done by introducing a perturbative 
expansion both in the asymmetry of the free energies and in the pulling speed. Moreover, 
our approach makes it possible to identify a critical rate that separates two well-dehned 
regimes. To the lowest order, this critical velocity has a linear dependence on the 
asymmetry of the potential. In spite of the crude approximations, our theory compares 
quite well with the numerical data even beyond the applicability regime. Moreover, 
our results provide a guide to interpret some inversions observed in the sequence of 
unfolding of the stable regions of the maltose-binding protein during its mechanical 


denaturation 17,23 


It must be stressed that to the lowest order in our theory, the system is sweeping 
the stationary branches of the force-extension curve. In this sense, the pulling process 
is very slow or adiabatic. Despite this adiabatic nature of the pulling process, it is not 
always the weakest unit to unfold hrst. As long as the pulling speed Up 7^ 0, the closer 
to the pulled terminal one unit is, the larger the force acting on it. This gradient in 
the distribution of the force across the protein, which increases with the pulling speed. 
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makes it possible that the last unit reaches first its limit of stability. 

We have limited ourselves to the investigation of the first unfolding event. However, 
our argument can be easily generalized to the next unfolding event: the difference is 
that the zero-th order approximation is no longer given by all the units sweeping the 
all-units-folded branch but by the sweeping of the branch with one module unfolded 
and the remainder folded. Then, a similar perturbative expansion around this zero-th 
order solution in powers of the asymmetry and the pulling speed would give the next 
unit that opens. 

If the biomolecule comprises several perfectly identical units, the asymmetry 
correction vanishes, because dUi = 0 (and thus dfi = 0) for all the units. In that 
case, our theory predicts that it is always the pulled unit that unfolds first. However, 
even in engineered modular proteins, slight differences from module to module may be 
present. In fact, this has lead to the analysis of the impact of quenched disorder in 
the force-extension curves of biomolecules 27 . In the present context, we may also 


introduce stochastic free energy deviations, following a certain probability distribution. 
Next, our theoretical approach can be applied to this system with quenched disorder 
in the free energies. Interestingly, evidence of dynamical disorder, that is, a fluctuating 
environment, has been recently brought to bear in stretching experiments 47 . The 


analysis of this situation needs a more complex theory, in which the free energy landscape 
fluctuates in time, and is outside of the scope of the present paper. 

There are methods that extract the free energy landscape from experimental data 
of pulling experiments, even when there are intermediates 36 40 . The resulting free 


energy is usually calculated as a function of the end-to-end distance of the molecule. 
Nevertheless, in order to apply our theory, we do not need this global energy landscape as 
a function of the end-to-distance of the molecule but each unit’s contribution thereto as a 
function of its own extension. In this regard, it is relevant to note that a similar velocity 
dependent unfolding pathway should also be found in modular proteins, although it 
has not been experimentally investigated to the best of our knowledge. In fact, we 
have analyzed a simple polyprotein model with a realistic potential, and observed a 
completely analogous behaviour. When all the modules are not identical, the weakest 
one will always open first for small enough pulling velocities. On the other hand, if the 
pulled unit is not the weakest, this will no longer be the case as the pulling speed is 
increased. Since the free energy of each module is experimentally accessible and the 
critical velocity lies on the experimental range, a reliable test of our theory could be 
done in modular proteins. Another possibility that deserves attention is to test our 
theory in ankyrin repeat proteins. In 48 , the unfolding of a consensus ankyrin repeat 
protein, NI6C, has been investigated. This protein is composed of eight repeats: the two 
capping repeats are different from the six identical central ones, that at the C-terminus 
(N-terminus) is weaker (stronger) than the rest. Pulling from the C-terminus (weakest 
unit), Lee et ah observe that the unfolding always starts from this end. This is what is 
expected from our theory, since when the pulled and the weakest units coincide, there is 
no competition between the asymmetry and kinetic terms. Therefore, it seems relevant 
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to carry out the same experiment but pulling from the N-terminus (strongest unit), in 
which a much richer phenomenology is to be expected. 

Very recently, sequential unfolding has been reported in a simple model 


35 


which makes it possible to understand the stepwise unfolding observed in force-clamp 
experiments with modular proteins 41-43 . This sequential unfolding appeared as 


a consequence of a depinning transition near the stability threshold introduced by 
the coupling between nearest neighbour units. Interestingly, the unfolding of the 
ankyrin repeat protein in does not only start from a well-defined end but it is 
also sequential, which may hint at the significance of this kind of short-ranged couplings 
in the experiment. Then, it seems also relevant to analyze whether a similar sequential 
unfolding is present in the model developed here, when a similar short ranged interaction 
between neighbouring units is considered. 
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Appendix A. Stability threshold 

To first order in the extension xt^b such that U'f{xi^b) = 0 verifies 

f^'"(4)(a;i,b-4)+e5/'(4) = 0, 

that is. 






(A.l) 


(A.2) 


The corresponding force at the stability threshold is obtained from (10). To the lowest 
order in the deviations. 


F,,b = F[{x,^b) ~ F[{lb) = Fb + ^SfM, 


(A.3) 


because the next term, U"'{ib){xi^b — is of the order of Therefore, the Tth 

module reaches its limit of stability at the time for which Xi = xf'^ -f fdxi = Xi^bi that 
is, when the length per unit ^ has the value verifying 


+ e 


5/(^i) - Wi 


or, equivalently. 


= e 


SMii 


5f{l 


U"{ii 



(A.4) 

^V'"(4)' 

(A.5) 


We know that £ ^ £b when ^ 0. But U'\£b) = 0 and thus we cannot substitute 

£b on the rhs of (A.5). On the other hand, this means that the dominant balance for 
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0 involves the Ihs and the hrst term on the rhs of (A.5). Therefore, making use of 
f/"(£)~f/'"(4)(^-4), weget 


(ii 




Smb) - 5f{h) 

U'"{h) ' 


(A.6) 


Since U'"{ib) < 0 (see hgure |^, this means that only the units with Sfi{ib) smaller 
than the average (that is, weaker than average) reach the limit of stability in the limit 
as Up —)■ 0. In fact, it is the weakest unit, that is, the unit with smallest 6fi{£b), that 
unfolds hrst. 


It is interesting to note that, in order to obtain (A.6), we have completely neglected 


the last term on the rhs of (A.5). Since, in turn, this term stems from the last term 


on the rhs of (A.2), to the lowest order we are solving the equation Xi = ib- In other 


words, to the lowest order the stability threshold can be considered to be given by the 
non-disordered, zero-asymmetry case, free energy U{x). For the sake of concreteness 
and simplicity, we have stuck to the asymmetry contribution dxi in this appendix, but 
the same condition Xi = ib would still be valid, had we taken into account the kinetic 


contribution Axi derived in section 3.2 The reason is that there is also a factor U"{i) 


in the denominator of Ax^, see (23a), and thus both the terms coming from 6xi and Ax* 


are dominant against the last term on the rhs of (A.2) 


Appendix B. Discrete inhomogeneous diffusion equation 

In this appendix, we briefly discuss a general procedure which is useful to solve linear 


difference equations similar to those in (19) and (22). The methods for solving difference 


equations often resemble those used for solving analogous differential equations; the 


latter may be thought of as the continuous limit of the former. Both (19) and (22) 


belong to the following general class of second-order linear difference equations for |/j, 

yi= 9 {y 2 ,---,yN-i), (B.ia) 

Vi+i + Vi-i -2yi = K, 1 <i < N, (B.lb) 

1/iv = h(?/2, • ■ ■ ,1/iv-i), (B.lc) 


in which g, h are arbitrary functions and A is a given constant. (B.lb) is a second-order 


linear difference equation, and (B.la) and (B.lc) are its boundary conditions. It may be 


thought of as a discrete inhomogeneous diffusion equation: j/j+i — j/j is the hrst discrete 
derivative, so that yi+i + yt-i — 2yi is the second discrete derivative [^. In complete 
analogy with the corresponding differential equation y" = K, the general solution of 


(B.lb) is 


• ^-2 

yi = co + cii + —I 


(B,2) 


in which cq and Ci are two arbitrary constants. The solution of (B.l) is, as usual. 


univocally determined by the boundary conditions, from which specihc values for cq 
and Cl are obtained. 
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Let us show that both (19) and (22) can be cast in the above form. Firstly, in (19b), 


it is easily seen that if we define j/j = 6xi + 6fi/U”{i), (B.lb) is obtained with K = 0. 
Secondly, in (22b), it is straightforward to identify yi = Axi and K = 


A simple calculation gives the constants cq and ci in (B.2) for each case, and thus the 
expressions for dxi and Axi in the main text. 


Appendix C. Order of the critical velocities 

Here, we prove that It is easy to show that this inequality follows if we 

have that 

xr Sfa,+AQi^k+2 - Vk) - - i^k+i) 

^Uk+2Pb) > - 


k'k+i — k'k 


(C.l) 


in which = 0^(0^ — 1). Due to (31), a^+i minimize Therefore, in particular. 


< ^"'“(0:^+2), which is readily shown to be equivalent to (C.l). 
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